% replication_code_table_4.m
% 
% This matlab code replicates table 4 in Miyamoto, W., T.L. Nguyen, and 
% H. Oh (2023), "In Search of Dominant Drivers of the Real Exchange Rate," 
% Review of Economics and Statistics.
% 


%% setup
clear all;close all;clc
setup_workspace;
setup_seed;

%% load dynare
% load data from dynare
cd(ps.fpath.Main);
cd(ps.fpath.Models);
load('gnk2_loop4_fevd_results.mat') 
load('gnk2_loop4_fevd_oo_.mat') 
cd(ps.fpath.Main);


%% compute FEVD of model's financial shock

% choose variables
var_name={'a1','a2','y1','y2','c1',... 
          'c2','h1','h2','I1','I2',...
          'nxy','qD','pi1','pi2','r1',...
          'r2','wedge','psi1','mo1','mo2',...
          'a1_a2','y1_y2','c1_c2','h1_h2','r1_r2',...
          'pi1_pi2','I1_I2'};

% choose shock
%%%var_shock=['e_a1']; % productivity
%%%var_shock=['e_mo1']; % monetary
var_shock=['e_psi1']; % financial

for i=1:length(oo_cell)
    eval(['fevd(:,:,:,i)=oo_cell{i}.conditional_variance_decomposition' ';'])
end

s_case=6; % case of the model

% FEVD: relative TFP (a1-a2)
var_pos = strmatch(var_name{21},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_TFP_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store TFP

% FEVD: relative output (y1-y2)
var_pos = strmatch(var_name{22},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_Y_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store Y

% FEVD: relative consumption (c1-c2)
var_pos = strmatch(var_name{23},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_C_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store C

% FEVD: relative hours worked (h1-h2)
var_pos = strmatch(var_name{24},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_H_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store H

% FEVD: net-export-to-output ratio (nxy)
var_pos = strmatch(var_name{11},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_NXY_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store NXY

% FEVD: real exchange rate (qD)
var_pos = strmatch(var_name{12},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_qD_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store qD

% FEVD: relative inflation (pi1-pi2)
var_pos = strmatch(var_name{26},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_PI_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store PI

% FEVD: relative nominal interest rate (r1-r2)
var_pos = strmatch(var_name{25},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_INT_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store interest rate

% FEVD: UIP wedge (wedge)
var_pos = strmatch(var_name{17},M_.endo_names,'exact');          
shk_pos = strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_wedge_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store wedge

% FEVD: relative investment (I1-I2)
var_pos=strmatch(var_name{27},M_.endo_names,'exact');          
shk_pos=strmatch(var_shock,M_.exo_names,'exact');
FEVDmodel_I_wedge = squeeze(fevd(var_pos,:,shk_pos,s_case)); % store wedge



%% load VAR estimation result using model's simulated data
load(['all_workspace'])
setup_workspace;

%% obtain FEVD from the VAR using model's simulated data
close all; clear tmp;
s_case=6; % case of the dynare model

% 1. FEVD of dominant relative output shock 
s_mbc=1; % target variable relative output=1

% get FEVD draws
for i=1:ps.nsim
    tmp(:,:,i)=Y_fevd_draws_all{s_case,i,s_mbc}; % case by simulation by mbc variable
end

% percentile
tmp50=prctile(tmp,50,3);
tmp10=prctile(tmp,10,3);
tmp90=prctile(tmp,90,3);
tmp16=prctile(tmp,16,3);
tmp84=prctile(tmp,84,3);


% FEVD: relative output
FEVDsim1_Y_50 = tmp50(1,:); FEVDsim1_Y_16 = tmp16(1,:); FEVDsim1_Y_84 = tmp84(1,:);
% FEVD: relative consumption
FEVDsim1_C_50 = tmp50(2,:); FEVDsim1_C_16 = tmp16(2,:); FEVDsim1_C_84 = tmp84(2,:);
% FEVD: relative hours worked
FEVDsim1_H_50 = tmp50(3,:); FEVDsim1_H_16 = tmp16(3,:); FEVDsim1_H_84 = tmp84(3,:);
% FEVD: NXY
FEVDsim1_NXY_50 = tmp50(4,:); FEVDsim1_NXY_16 = tmp16(4,:); FEVDsim1_NXY_84 = tmp84(4,:);
% FEVD: real exchange rate
FEVDsim1_qD_50 = tmp50(5,:); FEVDsim1_qD_16 = tmp16(5,:); FEVDsim1_qD_84 = tmp84(5,:);
% FEVD: relative inflation
FEVDsim1_PI_50 = tmp50(6,:); FEVDsim1_PI_16 = tmp16(6,:); FEVDsim1_PI_84 = tmp84(6,:);
% FEVD: relative nominal interest rate
FEVDsim1_INT_50 = tmp50(7,:); FEVDsim1_INT_16 = tmp16(7,:); FEVDsim1_INT_84 = tmp84(7,:);
% FEVD: UIP wedge (implied)
FEVDsim1_wedge_50 = tmp50(9,:); FEVDsim1_wedge_16 = tmp16(9,:); FEVDsim1_wedge_84 = tmp84(9,:);
% FEVD: relative investment
FEVDsim1_I_50 = tmp50(8,:); FEVDsim1_I_16 = tmp16(8,:); FEVDsim1_I_84 = tmp84(8,:);


% 2. FEVD of dominant RER shock 
clear tmp;
s_mbc=5; % target variable RER=5

% get FEVD draws
for i=1:ps.nsim
    tmp(:,:,i)=Y_fevd_draws_all{s_case,i,s_mbc}; % case by simulation by mbc variable    
end

% percentile
tmp50=prctile(tmp,50,3);
tmp10=prctile(tmp,10,3);
tmp90=prctile(tmp,90,3);
tmp16=prctile(tmp,16,3);
tmp84=prctile(tmp,84,3);


% FEVD: relative output
FEVDsim5_Y_50 = tmp50(1,:); FEVDsim5_Y_16 = tmp16(1,:); FEVDsim5_Y_84 = tmp84(1,:);
% FEVD: relative consumption
FEVDsim5_C_50 = tmp50(2,:); FEVDsim5_C_16 = tmp16(2,:); FEVDsim5_C_84 = tmp84(2,:);
% FEVD: relative hours worked
FEVDsim5_H_50 = tmp50(3,:); FEVDsim5_H_16 = tmp16(3,:); FEVDsim5_H_84 = tmp84(3,:);
% FEVD: NXY
FEVDsim5_NXY_50 = tmp50(4,:); FEVDsim5_NXY_16 = tmp16(4,:); FEVDsim5_NXY_84 = tmp84(4,:);
% FEVD: real exchange rate
FEVDsim5_qD_50 = tmp50(5,:); FEVDsim5_qD_16 = tmp16(5,:); FEVDsim5_qD_84 = tmp84(5,:);
% FEVD: relative inflation
FEVDsim5_PI_50 = tmp50(6,:); FEVDsim5_PI_16 = tmp16(6,:); FEVDsim5_PI_84 = tmp84(6,:);
% FEVD: relative nominal interest rate
FEVDsim5_INT_50 = tmp50(7,:); FEVDsim5_INT_16 = tmp16(7,:); FEVDsim5_INT_84 = tmp84(7,:);
% FEVD: UIP wedge (implied)
FEVDsim5_wedge_50 = tmp50(9,:); FEVDsim5_wedge_16 = tmp16(9,:); FEVDsim5_wedge_84 = tmp84(9,:);
% FEVD: relative investment
FEVDsim5_I_50 = tmp50(8,:); FEVDsim5_I_16 = tmp16(8,:); FEVDsim5_I_84 = tmp84(8,:);





%% summarize the FEVDs
s_time=0:20; % period for FEVD
T_fevd = length(s_time);

% FEVD of model's financial shock
num_modeluwedge = 100*[FEVDmodel_Y_wedge(1:T_fevd); ...
    FEVDmodel_C_wedge(1:T_fevd) ; ...
    FEVDmodel_I_wedge(1:T_fevd); ...
    FEVDmodel_H_wedge(1:T_fevd); ...
    FEVDmodel_NXY_wedge(1:T_fevd); ...
    FEVDmodel_qD_wedge(1:T_fevd); ...
    FEVDmodel_PI_wedge(1:T_fevd); ...
    FEVDmodel_INT_wedge(1:T_fevd)];


% FEVD of dominant relative output shock using simulated data

% 1. Median:
num1_simu50 = 100*[FEVDsim1_Y_50(1:T_fevd); ...
    FEVDsim1_C_50(1:T_fevd); ...
    FEVDsim1_I_50(1:T_fevd); ...
    FEVDsim1_H_50(1:T_fevd); ...
    FEVDsim1_NXY_50(1:T_fevd); ...
    FEVDsim1_qD_50(1:T_fevd); ...
    FEVDsim1_PI_50(1:T_fevd); ...
    FEVDsim1_INT_50(1:T_fevd)];
    
% 2. 16 percentile
num1_simu16 = 100*[FEVDsim1_Y_16(1:T_fevd); ...
    FEVDsim1_C_16(1:T_fevd); ...
    FEVDsim1_I_16(1:T_fevd); ...
    FEVDsim1_H_16(1:T_fevd); ...
    FEVDsim1_NXY_16(1:T_fevd); ...
    FEVDsim1_qD_16(1:T_fevd); ...
    FEVDsim1_PI_16(1:T_fevd); ...
    FEVDsim1_INT_16(1:T_fevd)];
    
% 3. 84 percentile
num1_simu84 = 100*[FEVDsim1_Y_84(1:T_fevd); ...
    FEVDsim1_C_84(1:T_fevd); ...
    FEVDsim1_I_84(1:T_fevd); ...
    FEVDsim1_H_84(1:T_fevd); ...
    FEVDsim1_NXY_84(1:T_fevd); ...
    FEVDsim1_qD_84(1:T_fevd); ...
    FEVDsim1_PI_84(1:T_fevd); ...
    FEVDsim1_INT_84(1:T_fevd)];


% FEVD of dominant RER shock using simulated data

% 1. Median:
num5_simu50 = 100*[FEVDsim5_Y_50(1:T_fevd); ...
    FEVDsim5_C_50(1:T_fevd); ...
    FEVDsim5_I_50(1:T_fevd); ...
    FEVDsim5_H_50(1:T_fevd); ...
    FEVDsim5_NXY_50(1:T_fevd); ...
    FEVDsim5_qD_50(1:T_fevd); ...
    FEVDsim5_PI_50(1:T_fevd); ...
    FEVDsim5_INT_50(1:T_fevd)];

% 2. 16 percentile
num5_simu16 = 100*[FEVDsim5_Y_16(1:T_fevd); ...
    FEVDsim5_C_16(1:T_fevd); ...
    FEVDsim5_I_16(1:T_fevd); ...
    FEVDsim5_H_16(1:T_fevd); ...
    FEVDsim5_NXY_16(1:T_fevd); ...
    FEVDsim5_qD_16(1:T_fevd); ...
    FEVDsim5_PI_16(1:T_fevd); ...
    FEVDsim5_INT_16(1:T_fevd)];
    
% 3. 84 percentile
num5_simu84 = 100*[FEVDsim5_Y_84(1:T_fevd); ...
    FEVDsim5_C_84(1:T_fevd); ...
    FEVDsim5_I_84(1:T_fevd); ...
    FEVDsim5_H_84(1:T_fevd); ...
    FEVDsim5_NXY_84(1:T_fevd); ...
    FEVDsim5_qD_84(1:T_fevd); ...
    FEVDsim5_PI_84(1:T_fevd); ...
    FEVDsim5_INT_84(1:T_fevd)];

    

%% display table 4 

% get FEVD for h=4 and h=20:
disp_num_modeluwedge = num_modeluwedge(:,[5 21])';

disp_num1_simu50 = num1_simu50(:,[5 21])';
disp_num1_simu16 = num1_simu16(:,[5 21])';
disp_num1_simu84 = num1_simu84(:,[5 21])';

disp_num5_simu50 = num5_simu50(:,[5 21])';
disp_num5_simu16 = num5_simu16(:,[5 21])';
disp_num5_simu84 = num5_simu84(:,[5 21])';


% create FEVD table
fevdtable_yshock_sim_h4    = array2table([disp_num1_simu50(1,:); disp_num1_simu16(1,:); disp_num1_simu84(1,:)], ...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate'},...
    'RowName',{'Median','16%','84%'});
fevdtable_yshock_sim_h20   = array2table([disp_num1_simu50(2,:); disp_num1_simu16(2,:); disp_num1_simu84(2,:)], ...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate'},...
    'RowName',{'Median','16%','84%'});

fevdtable_rershock_sim_h4  = array2table([disp_num5_simu50(1,:); disp_num5_simu16(1,:); disp_num5_simu84(1,:)], ...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate'},...
    'RowName',{'Median','16%','84%'});

fevdtable_rershock_sim_h20 = array2table([disp_num5_simu50(2,:); disp_num5_simu16(2,:); disp_num5_simu84(2,:)], ...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate'},...
    'RowName',{'Median','16%','84%'});

fevdtable_uipshock_model_h4 = array2table([disp_num_modeluwedge(1,:)], ...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate'},...
    'RowName',{'-'});

fevdtable_uipshock_model_h20 = array2table([disp_num_modeluwedge(2,:)], ...
    'VariableNames',{'Output','Consumption','Investment','Hours Worked','NXY','RER','Inflation','Interest Rate'},...
    'RowName',{'-'});

% display table 4
disp('-------')
disp('Table 4')
disp('-------')
disp('FEVD of dominant relative output shock; h=4')
disp(fevdtable_yshock_sim_h4)
disp(' ')
disp('FEVD of dominant relative output shock; h=20')
disp(fevdtable_yshock_sim_h20)
disp(' ')
disp('FEVD of dominant RER shock; h=4')
disp(fevdtable_rershock_sim_h4)
disp(' ')
disp('FEVD of dominant RER shock; h=20')
disp(fevdtable_rershock_sim_h20)
disp(' ')
disp('FEVD of model financial shock; h=4')
disp(fevdtable_uipshock_model_h4)
disp(' ')
disp('FEVD of model financial shock; h=20')
disp(fevdtable_uipshock_model_h20)
disp(' ')


